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It is shown how to formulate the ubiquitous quantum chemistry problem of calculating the ther- 
mal rate constant on a quantum computer. The resulting exact algorithm scales exponentially faster 
with the dimensionality of the system than all known "classical" algorithms for this problem. 
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I. INTRODUCTION 

It is well known that an exact calculation of the ther- 
mal rate constant is a problem that scales exponentially 
with the number of degrees_pf freedom.El But what if we 
had a quantum computerBEl (QC) at our disposal? As 
we will show here, then the calculation can be speeded 
up exponentially. Exponential speed ups on QCs have 
been demonstrated in a number of problems, the most 
famous of rwhich are Shor's algorithm for the factoring 
problempQ and Grover's algorithm for database search.0 
In the context of physics problems, exponential speed up 
has been demonstrated mostly in the context ojLsimula- 
tion of the many-body Schrodinger equationflu Other 
physics applications have_also been proposed, such as 
studying quantum chaosEl and Ising spin glasses.Eil At 
the experimental level QCs are still in a stage of in- 
fancy, although very impressive first steps towards irtu 
plementation have been taken using ions is, ion traps ,113 
atoms in high-finesse, microwave cavities,H3 and molec- 
ular spins in NMRo In particular, Chuang et al, us- 
ing a chloroform NMR-QC, recently implemented for the 
first time a quantum algorithm (Grover's) which outper- 
forms_any classical algorithm designed to solve the same 
task.Ej Another algorithm for which QCs offer an expo- 
nential speed up compared to classical computers, known 
as "Deutsch's problem, "£3 has also been implemented on 
a chloroform NMR. QCjlZl and by Jones et al. on a cyto- 
sine NMR-QC.Il§t2l These impressive achievements signal 
quite clearly, albeit for very simple applications at this 
point, that QCs may perhaps sooner than expected play 
an important role in simulations. For a comprehensive 
introduction to quantum computation, we refer the in- 
terested reader to a number of recent reviews .Ej 

The extensive body of work in quantum computation 
has, however, to date not addressed a computational 
problem of direct relevance to the quantum chemistry 
community (apart, of course, from the general simu- 
lation of the Schrodinger equation). In this work we 
show how to formulate the ubiquitous quantum chem- 
istry problem of calculating the thermal rate constant 
on a QC, and by doing so, how the calculation can be 
speeded up exponentially with the number of degrees of 
freedom. The rate constant is the single most important 
number characterizing chemical reactions, and thus great 
efforts have been invested in designing efficient and exact 



"classical" computational ways to obtain it. Important 
progress along this line has been made by Miller et al., Ell 
Light et a/.,E3, and Manthe et al.jEi based on|-tJa6j effi- 
cient evaluation of the flux correlation function.EJE3 Ap- 
proximate methods have also been developed for obtain- 
ing the rate constants. Far example, the popular mixed 
quantum-classical modelE£rE3 whereby one integrates the 
time-dependent Schrodinger equation for (a few) degrees 
of freedom that are treated quantum mechanically, simul- 
taneously with the classical equations of motion for the 
(many) degrees of freedom that are treated by classical 
mechanics; the semiclassical initial value repsesentationEa 
(SC-IVR) that has had a rebirth of interestE3~Ell as a way 
for including quantum effects in molecular dynamics sim- 
ulations; and the (further) linearizing approximation to 
the SC-IVR which leads to a much simpler form for the 
rate expression£3 While these methods enjoy favorable 
computational scaling properties, they are inherently ap- 
proximate and thus not in the context of the present 
paper. 

In spite of these significant advances, exact classical 
algorithms can at most achieve a polynomial speed up 
in a problem that inherently scales exponentially. In- 
deed, when classical algorithms are described as 0(N 3 ) 
instead of 0(N 2 ), it should be remembered that N is 
itself exponentially large! The anticipated advance in 
quantum computation would therefore have revolution- 
ary consequences for quantum chemistry, rendering "clas- 
sical" simulation methods essentially obsolete. Of course, 
QCs are still at best many years away from reaching the 
point of replacing classical computers. Nevertheless, it 
is of considerable interest to exhibit an explicit QC al- 
gorithm for a problem as central as the computation of 
the rate constant, and this is the task we undertake here. 
The paper is organized as follows: In Sec. || we intro- 
duce some pertinent concepts of quantum computation. 
Then, in Sec. [II we briefly rederive the exact quantum 
expression for the rate constant. The next sections are 
the heart of the paper, where the QC algorithm for calcu- 
lation of the rate constant is described in detail. Sec. VI 
concludes. 
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II. A BRIEF INTRODUCTION TO QUANTUM 
COMPUTATION 



Let us denote all the degrees of freedom of the prob- 
lem at hand by the collective variable q = {q\, qja)> 
and let us assume for simplicity that after discretization 
the number of points per degree of freedom is 2 l . Then 
the Hilbert space for a given Hamiltonian H (q, p) is of 
dimension 



V 



V 



N = 2 V , 



IM 



(1) 



whence the exponential scaling. To set up the prob- 
lem on a QC one introduces a "register" of v "qubits" 
(two- level systems), which can be in a superposition state 
\4>i) = ai\0i) +bi\li) (with H 2 + |^| 2 = 1). Each group of 
I qubits corresponds to one of the degrees of freedom qj . 
Initially the quantum register |$) is in a direct product 
state with all qubits in the |0) state: 



1$) 



(2) 



Allowed operations on the register are all the unitary 
transformations (corresponding to propagation of the 
register), and all measurements, i.e., projections onto 
subspaces of the full register Hilbert space Tt. However, 
by convention the unitary transformations should be ex- 
plicitly given in terms of operations on single and two 
qubits at the most, since it is such one and two-qubit 
"gates" that one can expect to construct in practice. 
Also, allowing arbitrarily large gates would not consti- 
tute a general-purpose computer (this is similar to the 
situation with classical computers, where the number of 
distinct logical elements is a small and finite set). Sim- 
plify ine_an early construction by Deutsch,cl it has been 
provenEj that the set of all single-qubit gates [the group 
U(2)} with in addition the " controlled-nof (CNOT) 



|ei,e a ) ' — ► |ei,(ei+e 2 )mod2) (e, = 0, 1) , 



(3) 



form a universal set of gates: all unitary transformations, 
of arbitrary size, can be constructed using a polynomi- 
ally large set of the 1- and 2-qubit gates. Therefore the 
restriction to this set of gates is sufficient to simulate 
any computable function. This, of course, includes (and 
potentially exceeds) anything that is computable on a 
classical computer. 

Now, let us see how to set up a superposition state on 
the quantum register corresponding to all possible initial 
classical positions. In the basis |0) = (J) and |1) = (J), 
one applies the one-qubit unitary "Hadamard transform" 



(4) 



which is a ir/2 rotation on the Bloch sphere of each qubit. 
Thus: 



|$>^|*')=® R l i>=QS> 

t=l i=l 
JV-1 



V2 



(5) 



3=0 



where j is the decimal representation of the register state 
with the corresponding binary value, so |j) is a conve- 
nient shorthand notation for a tensor product of v single- 
particle states. Eq.(||) represents the desired superposi- 
tion over all initial positions. For example, let I = 2 (so 
the number of grid points per degree of freedom is 4) 
and M = 2 (e.g., q\ = x and q% = y va. & 2-dimensional 
problem involving a linear triatomic vibrating molecule) 
so N = 16; how is q — (2Ax, 3Ay) represented in the 
register? In binary, q — ({1,0} Ax, {1,1} Ay), so that 
j = 11, i.e., \j) = |1) ® |0) ® |1) ® |1), corresponds to 
x = 2Ax and y = 3Ay. In this way the quantum reg- 
ister supports a superposition over all discretized values 
of the degrees of freedom q. By linearity, unitary evolu- 
tion of the register amounts to parallel propagation on 
all of the exponentially many grid points. We will find it 
convenient to work from now on with j as the collective 
degrees of freedom variable, instead of q. The transfor- 
mation between the two is straightforward. Note further 
that the register states \j) are position eigenstates. 

Suppose now that the initial wavefunction |^(0)) in 
our scattering problem has amplitude cxj to be at po- 
sition j. As explained below, it is in fact not essential 
to utilize this initial condition, and an equal superposi- 
tion of all possible position eigenstates will be sufficient 
in most cases. Nevertheless, as shown by Zalkatl (and 
will not be repeated here), it is possible to initialize the 
register to the state 



1*0 



JV-l 
j=0 



(G) 



by means of a suitable unitary transformation. This then 
represents |^(0)) on the QC. The dynamics in the scat- 
tering problem is determined by the unitary propagator 
U = e - lHt / h : = U|*(0)>. The crucial advantage 

offered by a QC is that, as will be shown below, it is 
possible to efficiently implement this propagator on the 
quantum register, so that: 



I*") 



(7) 



In this way we have set up a 1-1 correspondence between 
the QC (| < I > )) and the dynamics of the problem of interest 
(|\1/)). All the relevant information can be extracted from 
this simulation by observing the states of the qubits. 
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III. THE THERMAL RATE CONSTANT VIA THE 
FLUX CORRELATION FUNCTION FORMALISM 

Let us now turn to the scattering problem and define 
a flux operator 



F=-[H,h[ S (g)]] , 



(8) 



where h is the Heaviside function and the condition 
s(q) = defines the dividing surface. 

The thermal rate constant is written as the-time inte- 
gral of the flux-flux autocorrelation functions 



k(T) 



1 



Qr(T) J 



dtC f (t) , 



(9) 



where Q r (T) is the reactant partition function per unit 
volume, and: 

C f (t) = Irle-WFe-We^Fe-™^] 

_ Tr jp e iHt/?i-/3H/2 p e -iHt/ft-/3H/2j 

= Tr[Fe iHT */ R Fe- ,Hr / fi ] (10) 

and where we have also defined for convenience the com- 
plex "time" 



T = t-ih/3/2. 



(11) 



Evaluating the trace in the energy eigen-basis {\n}} (with 
H|n) = E n \n)) we obtain: 

C f (t) =^(n|Fe iHT *- /n Fe- iHr / R |n) 

n 

= ^(n\Fe tUT */ h |m)(m|Fe- iHT / r » 

= J2 e iE m T*/H e -iE n r/h ( n | F | m )( m | F | n ) 

n,m 

= J2 e -^ +E -^ 2 e l ^- E ^ h \{n\F\m)\ 2 . (12) 
Using the commutator form for F, Eq.(||), we find: 



(n\F\m) = ~<n| [Hh[s(q)} - h[ s (q)}H] \m) 
= UE n -E m )(n\h[s(q)]\m). 



(13) 



Recall from Eq.(]|) that the quantum register naturally 
supports a superposition over position eigenstates. Ac- 
cordingly, let us represent \n) in the discretized position 
basis which we from now on identify with the QC's 

"computational basis" {\j}} (indeed, the correspondence 
is 1-1). Thus, let us expand the energy eigenstates as: 



JV-l 



(14) 



Clearly, the Heaviside function h [s(q)] is diagonal in this 
basis, so that from Eq. ([l3|) : 



N-l 



(n\F\m) =jr(E n - E m ) ^ a*{n) aj (m)h [s(j)] . (15) 



Hence, finally: 

cm = £ E 

N-l 



e -P(E m +E n )/2 e i(E m -E n )t/K ^ _ Em) 2 



aj(n)aj(m)h[s(j)] 

3=0 



(16) 



Our task is therefore to find an algorithm that calcu- 
lates the spectrum {E n }, and the position amplitudes 
{Oj(n)}^) 1 for each of the eigenstates \ j). With these in 



hand the rest of the calculation [summations in Eq.(|l6|)l 
can be efficiently implemented on a classical computer. 

At first sight it appears that there are two problems: 
(1) The spectrum may contain exponentially many ener- 
gies; (2) The number of measurements needed must be 
exponentially large (even if the number of energy eigen- 
states is polynomial) since there are exponentially many 
position eigenstates Regarding (1), the summation 
should indeed extend over a polynomially large number 
of energies only. This, however, is reasonable if the spec- 
trum is sufficiently degenerate and, more importantly, 
since the exponential decrease due to the Boltzman fac- 
tors will effectively eliminate the higher end of the spec- 
trum E > fcB? 1 +(barrier height) even if it is not degen- 
erate. As for (2), the question is how to get a reasonable 
sample of the distribution. For example, if the distribu- 
tion is badly behaved in position space, one could Fourier 
transform to momentum space and sample there. At any 
rate, this problem is identical to the one faced by any 
classical simulation, where one settles for a Monte Carlo 
sampling of the wavefunction. This is not the source 
of the classical bottleneck, and we will adopt the same 
Monte Carlo approach in the quantum simulation. Note, 
however, that unlike the classical case where attention 
has to be given to the issue of generating statistically 
independent Monte Carlo samples, the quantum simula- 
tion automatically generates_jtruly independent samples 
by the projection postulate.EJ 

How can we efficiently calculate the eigenstates and 
the spectrum? Solving the time- independent Schrodinger 
equation (SE) 



(17) 



can be done on a QC by transforming the problem to the 
time-dependent SE and propagating the dynamics with 
the unitary time-evolution operator U = e~ iH *. Each 
energy eigenvalue and eigenstate can then be obtained 
via known quantum algorithms, to be detailed below, in 
polynomial time. 
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IV. GENERAL OUTLINE OF THE ALGORITHM 

We now give the algorithm in general terms, to be de- 
fined more precisely in the next section. 

1. Prepare a register as in Eq.(||), and attach some 
"ancilla" qubits to it, also in the |0) state. These 
will serve as a quantum scratch-pad to record the 
results of intermediate measurements. From now 
on we will distinguish between the "main" and an- 
cillary registers. 

2. If a good guess for the initial wavefunction is 
known, initialize the register to it as in Eq.(|^). 
Else initialize the register to an equal superposi- 
tion. Since the computational basis states are po- 
sition eigenstates in all likelihood they are not en- 
ergy eigenstates, so will not be stationary under 
the SE dynamics. Thus except if the equal super- 
position corresponds to some undesirable position 
- such as very high above the barrier so that disso- 
ciation sets in immediately - it is as good a guess as 
any. In fact, any random (but reproducible) initial 
distribution will do. 

3. "Propagate" the register in parallel for a time t. 
This corresponds to a parallel evolution of all the 
position eigenstates. The propagation is done ve*j 
much in analogy to the classical FFT method. — 
in particular the split time propagation schemeE! 
Namely, the potential part is diagonal and can be 
implemented directly, whereas for the kinetic part 
it is necessary to Fourier transform to and back 
from momentum space. 

4. Perform a "von Neumann" measurement (see 
Sec. |ypj ) on the ancillary register using the Hamil- 
tonian (energy) as the observable. This accom- 
plishes a double purpose: 

(a) It allows to obtain an energy E n by measuring 
the ancillas. 

(b) It provides a means to sample the energy- 
position amplitudes aj(n). 

5. Repeat steps 1-4 many times until the distribution 
is converged to the desired accuracy for all relevant 
eigenstates. The number of required repetitions is 
proportional to this accuracy. 

6. Calculate (classically) the sums in Eq.([lG|). 



physics of this initialization step depends on the QC 
implementation. One conceivable way is cooling to the 
ground state. 



B. Inputting the Initial Wavefunction 

If necessary one inputs the initial wavefunction by the 
technique of Zalkafl Else one employs the Hadamard ro- 
tations technique to create an equal superposition over 
position states, as in Eq.(|J). In the former case the reg- 
ister will be in the state: 



' N-l 



\0i) 



(18) 



In the latter case all ctj = 1. 



C Quantum Propagation Algorithm 

This subsection is the heart of the algorithm; it builds 
on the approach of ZalkaJj Assume for simplicity that we 
have a single particle of mass m in an external potential 
V(q). The full Green's function for arbitrary time t is: 



G(xi,x 2 ;t) 



( Xl \e- iKt ' h \x 2 ) 



(19) 



For short time steps At <C 1 /E (E is a typical energy of 
the system) this becomes approximately: 



G(xi , x%\ At) = k exp 



(xi - x 2 y 

2At 



zV(xi)Ai 



(20) 



where k is a normalization factor. Applying this to the 
amplitudes is equivalent to acting on the basis states 
with the inverse transformation. Thus the position eigen- 
states, properly normalized, transform as: 



li) — u|j) = 

JV-l 

vw E ex p 

j'=a 



-im- 



(j-jfAx 2 



2At 



iV(jAx)At 



\f)- 
(21) 



This is carried out in parallel on the entire superposition 
SjLcj 1 li)- Suppose the time-step and spatial resolution 
are adjusted so that: 



V. THE ALGORITHM IN DETAIL 



A. Initialization 



Here the register is initialized to the state |<E>) = 
|0i)> where the last v qubits are ancillas. The 



At 



N 



(22) 



Then by expanding the exponent Eq.(pT|) can be writ- 
ten as a succession of a diagonal transformation, Fourier 
transform, and another diagonal transformation, all uni- 
tary: 
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U|i) = exp [iF 2 (j)] F(j,f) exp [tfitf')] \f) 
where: 



(23) 



Fl(j) = -7T 



r 

iV 

•2 



Fa(j) = -7r^ + VC7'Aa:)A* 



N-l 

E 



.7'=0 



exp 



2ttz 



I/) 



(24) 



Eq. ( p2| ) tells us how many qubits f = log 2 N are needed 
for given Ax and At: 



log 2 



2?r At 
to Ax 2 



(25) 



The special form of Eq. ( p3] ) , involving diagonal transfor- 
mations and a Fourier transform, is due to the structure 
of the Hamiltonian operator as a sum of operators diag- 
onal in coordinate and momentum space. As mentioned 
above, this is very similar to the situation that arises in 
the classical FFT method for solving the SE.EJ 



x = J2t=-k x $-> we have: |x) = |x_ fc ) ® \x-k+i) 
\xk-i), where xi — 0, 1. In the standard basis |0) = ( ) 
|1) = (,), consider the following unitary operation: 



Q 



fe-i 



l=-k 



1 

e l2 ' 



(26) 



The I th 2 x 2 matrix is a unitary operation in the Hilbert 
space of qubit number I. Thus: 



I p i )\xi) = e^\x l} 



Therefore the full result is: 
fe-i 



Q\x) 



e lXl2 \xi) = e 



Ek—1 
x l 
l = ~k 1 



k-l 



(27) 



\xi) = e ix \x) 



l=-k 



l=-k 



(28) 



as required. 



2. The Quantum Fourier Transform 



1. Diagonal Transformations: 

Consider first executing the diagonal unitary transfor- 
mations \j) i — > exp [iF(j)] \ which can be done as fol- 
lows, using the ancillary register |Eq.(|l8|)l, in the state 
|0) = (z$)j!L v |0i). The number v of qubits in this register 
depends on the accuracy with which F needs to be eval- 
uated (see immediately below) . Then the following steps 
are applied: 

1. \j, 0) i — ► \j,F(j)): evaluation of F and storage of 
the result in the ancillary register; 

2. \j, F(j)} i — ► exp [iF(j)} \j, F(j)) : introducing the 
phase; 

3. exp[iF(j)] \j,F(j)) i — > exp[iF(j)] \ j,0): inversion 
of step 1 in order to clear the ancillary register. 

Step 1 requires that it is possible to evaluate an arbi- 
trary function and store the result. This is very similar 
to the equivalent classical problem, for which algorithms 
are known using just the elementary classical gates. The 
same can be done in the quantum case, by breaking up 
the evaluation into elementary arithmetic operations, for 
which quantum algorithms have been designedlltH We 
will not dwell on this issue here. Step 3 is just the reverse 
of step 1 and can therefore be implemented by running 
the inverse unitary transformation. 

Step 2 has no classical analogue since it involves 
phases. It can be implemented if one knows how to do 
\x) i — ► exp [ix] \x). This can be done by simple single- 
qubit phase-shifts. Let v = 2k. Using a binary expansion 



The quantum Fourier transform, (QFT) algorithm has 
been discussed extensively,BH~c2l and some|— beautiful 
connections to group theory have been made.c2l In view 
of its central importance in the algorithm for solving the 
SE (and indeed in all efficient quantum algorithm found 
so far!), we present aJarief derivation here, using the ap- 
proach of Cleve et alzB 

The QFT was defined in Eq.(|||). Using the binary- 
decimal notation j/N = Q-jxfe—ju (recall that N = 2") 
where j± = 0, 1 etc., we note first that: 



o 2Tti(0.j 1 j 2 ...j^)jl\A' 



(29) 



It follows that: 



N-l 
£ eXP 

j'=a 



2m- 

N 



(|0) + e 2 «(°*-i>)|l)^ g ... 
® (|0)+e 27ri ( ^' 2 - i -)|l)) , (30) 



by expanding out the product on the right-hand-side 
and a term-by-term comparison. Thus the Fourier- 
transformed state in Eq.([30]) is in fact an "unentangled" 
direct product. This fact greatly simplifies the imple- 
mentation of the QFT. 

To perform the QFT, one first applies a Hadamard ro- 
tation [Eq.(||)] to \j\) (the first qubit of with the 
result: 



G 
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R|ji) = (|Q) + (-l)*|l» = (|Q) + e 2 ™^)|l>) , (31) 

so: \j) .— > (|0)+e 2 "(°^)|l))|j 2 ,...,^). Let us now 
define a new single-qubit operation, similar to Q from 
Eq.©: 



Q 



i — 



1 o 

e 2 "/ 2 ' 



(32) 



This operation is applied on the first qubit subject 
to a control by a second qubit \ji) (which itself does not 
change) : a "controlled rotation" . Namely, if ji = one 
does nothing, if it is 1, one applies Q;. This can be 
written as the following unitary transformation in the 4- 
dimensional Hilbert space of the two qubits, in the stan- 
dard basis = |00) = (1,0,0,0), |01) = (0,1,0,0), 

|10) = (0,0,1,0), |11) = (0,0,0,1): 



(33) 




CQi = 

After applying CQ2 one obtains: 

\ + e 2Ti(0.jij a )M\ 

Next a "controlled-Q3" is applied, yielding: 



10) 



(34) 



(35) 



Clearly, this process will eventually generate the desired 
phase in the superposition state of the first qubit [corre- 
sponding to the last qubit in Eq.(BG)]: 



U cc *n R 



|ji)= (|0)+e 27ri (°-^-^|l)) . (36) 



where the terms in the product from here onwards are 
applied tow index first. 

Now we turn to the second qubit. Again, a 
Hadamard rotation on it has the effect of: R|j 2 ) = 
(|0) +e 2 * li -°-^\l)). This is followed by a controlled- 
Q 2 , conditioned upon |j 3 ): (|0) + e 27r ^°*)|l)) 1 — ► 
(|0) + e 27ri (°- J ' 2j ' 3 )|l)). After the full operation on \ j 2 ) one 
obtains: 



n r 



j 2 



|i 2 ) = (|0)+e 2 ^°-^->)|l)) |j 2 >. 

(37) 



which corresponds to the one before last qubit in Eq. ( J3C|) . 

The method to generate the entire product in Eq.(|30|) 
should now be clear; collecting all the transformations 
yields: 



n 

P =i 



' v—p 



n CC M R 



S Z=2 



lii, -,ju) 



(|0) + ^(o-j'im-^)]!)) ® ... 

<8> (JO) +e 2 *- 1 ^|l)) ® (|0) 



O 27ri(0.>)| 



(38) 



Up to an unimportant bit reversal (which can easily be 
rectified by permuting the role of the qubits in the trans- 
formations above), this is exactly the desired result. In 
other words, the QFT is simply: 



0=1 



JJ CQi ] R 



vi=2 



(39) 



This will be applied in parallel, by virtue of the super- 
position principle, on all position eigenstates \j). Most 
importantly, the number of operations (single- and two- 
qubit) needed to implement the QFT is seen to be a mere 
v(y — l)/2. This is to be compared to the v2 u operations 
required classically, and as emphasized above, is the "se- 
cret" behind the quantum speedup. 



D. von Neuman Measurements 



Combining Eqs. 
is in the state 



and (|2l|), at this point the register 



i$"') = E^ u ii) = E^'(*)ij") 
3 3' 

i Jj ,(t)='£a j G- 1 (j,j , ;t). (40) 



A parallel propagation has occurred on all the position 
eigenstates. By measuring the qubits one by one, i.e., 
projecting onto a random position eigenstate and 
repeating this process many times while collecting the 
statistics, one can sample the electronic density function 
\ipj>(t)\ 2 . Our goal was to find the energ?/-spectrum and 
energy-position amplitudes a 3 - (n) , so these should be ob- 
tained from the simulation. This can be do«e using the 
so-called "von Neuman measurement" trick.cl We will re- 
quire an additional propagation step. 

A "measurement apparatus" that can be made to in- 
teract with the QC is introduced, and is assumed to be 
equivalent to a 1-dimensional quantum mechanical par- 
ticle. That is, its Hilbert space is spanned by the basis 
vectors |x), x real, with X|x) — x\x). In practice this will 
be another ancillary quantum register, consisting of, say, 
K qubits. Now, let us expand the position eigenstates 
in terms of the complete set of energy eigenstates 
[recall Eq.Q]: 



2^ a 3 



(n)|n) 



(41) 
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Consider next the joint evolution of an energy eigenstate 
\n) and the apparatus state \x) {x is arbitrary), under the 
unitary operator U = cxp(zH Pt/fi), where [X, P] = ifi. 
Here H acts on the main register and X, P act on the ap- 
paratus, so [X, H] = [P,H] = 0. We will shortly discuss 
the implementation of U. Consider first a formal Taylor 
expansion of exp(iH Pt/h), which yields: 



V\n)\x) 



oo 1 



1=0 



\n)\x + tE n 



(42) 



Thus U does not change the energy eigenstate, but has 
the effect of "shifting the dial x" by an amount propor- 
tional to the energy E n . The effect on the position eigen- 
state |/) will be: 

\J\ j')\x)=J2 4 (n)\n)\x + tE n ), (43) 

n 

and the effect on the full superposition of Eq.(|40|) is: 
U|$'") \x) = £ (t) E 4 (n)\n)\x + tE n ) 

j< n 

= Y / ^(t)\n)\x + tE n ) 

n 

= Y, a *A n )M*)- ( 44 ) 



Now suppose we observe the state of the apparatus. From 
Eq. (|i4|) it is clear that the apparatus has become entan- 
gled with the QC, and by performing the observation 
the superposition will collapse onto a particular state 
\m)\x + tE m ). This happens with probability |£ m (£)| 2 - 
Recall that \x + tE m ) is represented in binary by the 
qubits of the apparatus. Since t is a parameter of the 
simulation and x is known, all that remains is to mea- 
sure the apparatus qubit by qubit, to obtain the energy 
eigenvalue E m \ The accuracy with which these numbers 
are obtained is proportional to the number of simulation 
steps.Q 

To implement U it is necessary to Fourier transform 
the \x) register, just like in the classical FFT case. Specif- 
ically, let us define the Fourier transform pair: 



2 K -1 

T ' = ^ = T7w E e ~ lxv/U 



x) 



2^-1 



\p)- 



(45) 



p=0 



Then starting from the initial apparatus state \x), U can 
be implemented as follows: 



\n)\x) 



U 



\n) \p) 



i — ► In) 



p=0 

\n)\x + E n t), 



(46) 



in agreement with Eq. (|42j) . 



E. Extracting the Amplitudes from the 
Measurements 

Note further that after observation of the apparatus, 
the state of the main register has been projected onto 
|m), an energy eigenstate. The U propagation had a re- 
markable outcome: it transformed the information in the 
main register from a superposition over position eigen- 
states \j') to one over energy eigenstates \n). Unfortu- 
nately, apart from utilizing this in a subsequent evolu- 
tion, this, however, does not appear to be particularly 
useful, for it is not clear how the energy eigenstates are 
enumerated in the main register after the von Neuman 
propagation. Nevertheless, it is possible to obtain the 
amplitudes a,j{n) needed to complete the calculation in 
Eq.Q. Note first that by Eq.@: 



^fn(*)a,-(«) = 1pj(t). 



(47) 



Now, the simulation yields, by performing the whole 
procedure a sufficient number of times, an estimate of 
the probabilities \ipj{t)\ 2 [Eq.©] and |£™(i)| 2 [Eq.(§J)]. 
Thus to fully specify the complex numbers cij(n), it is 
necessary to also know their phases, as well as those of 
the £, n (t) and tpj{t). 

To obtain the phases, we note first that it is sufficient 
to know only the signs, since no generality is lost by 
employing a real initial wavefunction ^(O)). The signs 
can then be obtained with the help of a simple trick, 
which we will illustrate on a generic 2-qubit register state 
\ip) = a \00) + ai|01) + a 2 |10) + a 3 |ll). Given repeated 
preparations of this we perform the following set of 
measurements: 

• Observation of the two qubits in \ip). 

• A Hadamard transform on the first qubit, followed 
by observation of the two qubits. 

• A Hadamard transform on the second qubit, fol- 
lowed by observation of the two qubits. 



e lE " pt/h \n)\p) 



The first step yields an estimate 
one yields an estimate of |ao ± a\ \ 
der the Hadamard transform l^) 

(a -a 1 )\01) + (a 2 + a 3 )\W) + (a 2 - 
third step yields an estimate of 
since \ip) i — ► ^[(a + a 2 )|00) 4 

a 3 )|10) + (ai-a 3 )|H)]. Clearly, 
information for extraction of the 



of the | Oi | . The second 
and | a 2 ± 03 1 , since un- 
^^[(ao + «i)|00) + 
-03)|11)]. Similarly, the 
a ± a 2 | and |ai ± a 3 |, 
- Oo - a 2 )|01) + (ai + 

this provides sufficient 
signs of all amplitudes. 
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The generalization to a v-bit register is obvious: one per- 
forms Hadamard rotations on all v qubits. This then 
yields {|ao ± a±\, |o2 ± az\, \a>A ± asl, •■■} (after Hadamard 
on first qubit); {|ao ± aal, |ai ± 03!, |oi4 ± ae|, ■■■} (after 
Hadamard on second qubit); etc. After each Hadamard 
rotation there are 2 V coefficients to be estimated. This 
exponential "Monte-Carlo scaling" is the same as the one 
we encountered before and is not considered a slow-down 
for the reasons detailed above. The additional compu- 
tational cost is in the Hadamard rotations, v of which 
must be performed. This does therefore not affect the 
efficiency of the algorithm. At the end of the process, if 
the whole phase space has been sampled, one is left with 
v^y absolute values equations, which contain sufficient 
information to solve for the signs of all the amplitudes. 
In practice one will of course sample only a small (poly- 
nomial) portion of the phase space, and care must then 
be taken to obtain sufficient equations of the type above 
to uniquely determine the signs of the amplitudes of in- 
terest. 



F. Repetition 

The steps outlined above generate the energies {E n } 
and estimates of amplitudes {aj(n)} needed to perform 
the sum in Eq.@. The whole process must now be re- 
peated many times, on the order of the required accuracy, 
in order to complete the calculation. Due to the speed 
up in the implementation of the propagation step, the 
algorithm performs exponentially faster than any exact 
classical algorithm designed to solve the same task. 



VI. DISCUSSION AND CONCLUSIONS 

QCs are still far from being a panacea, and doubts 
have been raised whether they will ever replace ordinary, 
classical computers. Such worries are invariably based 
on the immense difficulties associated with maintain- 
ing phase coherence throughout the computation, i.e., 
the "decoherence problem." However, a remarkable the- 
ory of quantum error correction codes has recently been 
constructed, CJ in which a "logical qubit" is encoded in the 
larger Hilbert space of several physical qubits.L3 It has 
been shown that as long as the error rate is sufficiently 
small, it is possible to perform fault-tolerant quantum 
computation, i.e., the computation can be stabilized and 
be made fully robust to errors. o These advances greatly 
enhance the prospects of the eventual construction of use- 
ful QCs, beyond the current highly rudimentary proto- 
types. Building on these hopes, we have presented here 
an algorithm for calculating the thermal rate constant 
on a QC. The algorithm involves an initialization step 
of the QC into an equal superposition of position eigen- 
states; a propagation using an adaptation to QCs of the 
well-known FFT technique; and finally, a sequence of 



measurements yielding the energy spectrum and ampli- 
tudes. Under reasonable assumptions about the distri- 
bution of energy eigenvalues the algorithm runs in poly- 
nomial time. The algorithm thus outperforms any exact 
classical simulation, which is bound to be exponential. 
This clearly demonstrates the potential utility of QCs in 
future applications to quantum chemistry problems. 

Our approach was somewhat of a "brute force" one, 
in that we did not attempt to optimize the algorithm 
using such fruitful concepts as "direct and correct" low- 
rank expressions for the rate constant £1 Such optimiza- 
tions, while ineffectual in altering the essential exponen- 
tial speed up achieved by use of a QC, may still be im- 
portant in practice, especially in the early stages of the 
application on a small-scale QC of an algorithm such as 
described here. Further work is hence desirable to opti- 
mize the algorithm. 

Finally, it would be interesting to check the effect of 
noise and other types of errors affecting the evolution of 
the QC on the present algorithm. It has been shown, e.g., 
in the case of the ion trap QC, that factoring becomes 
impossible once random phase fluctuations in the laser 
pulses exceed a certain thresholdHa We intend to study 
similar noise-related issues using numerical simulations 
in the context of the present algorithm. 
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